Baltimore Ceasefire 365 is a city-wide call asking Baltimore residents to avoid having any murders through quarterly Ceasefires and Peace Challenges (Feb, May, August, and Nov). Here we model the impact of the ceasefire on shootings.

Shootings in Baltimore

Baltimore releases detailed data on issues relevant to the city, including crime. This allows us to get a pretty precise idea of the distribution of shootings in Baltimore.

Shootings per day seem to have increased somewhat over time, but the time series is complicated. 

library(tidyverse)
library(scales)

# data backed up on github
bpd <- read_csv("https://raw.githubusercontent.com/peterphalen/ceasfire/master/BPD_Part_1_Victim_Based_Crime_Data.csv")

bpd <- subset(bpd, Description == "SHOOTING")

bpd$CrimeDate <- as.Date(bpd$CrimeDate, format = "%m/%d/%Y")

# there are many crimes per day. collapse to daily counts
daily <- bpd %>% group_by(CrimeDate) %>% summarise(shootings = n())

# fill missing dates, because some had no shotings
full.ts <- data.frame(CrimeDate = seq(daily$CrimeDate[1], 
                                      daily$CrimeDate[nrow(daily)], by="day"))
daily <- full_join(full.ts,daily)
daily <- daily %>% group_by(CrimeDate) %>% mutate_all(funs(ifelse(is.na(.),0,.)))

ggplot(daily) +
  aes(x=CrimeDate, y=shootings) +
  geom_point(alpha=.2) + 
  xlab("date") +
  ylab("shootings") +
  scale_x_date(labels = date_format("%b %Y")) +
  ggtitle(" ", 
          subtitle="Baltimore (2012-present)")

These shootings disproportionately affect black communities. 

library(geojsonio)
library(leaflet)

bpd$Neighborhood <- factor(bpd$Neighborhood)
bpd <- subset(bpd, !is.na(Neighborhood))

count <- bpd %>%
  group_by(Neighborhood) %>%
  summarise(total.count=n()) 

nbds <- geojsonio::geojson_read("https://gis-baltimore.opendata.arcgis.com/datasets/1ca93e68f11541d4b59a63243725c4b7_0.geojson", what = "sp")

getcount <- function(neighborhood){
  nbd <- as.character(neighborhood)
  nbds <- as.character(count$Neighborhood)
  if(nbd %in% nbds){
    count <- count[count$Neighborhood == nbd,]$total.count
    return(count)
  }
  if(!(nbd %in% nbds)){
    return(0)
  }
}

nbds$count <- sapply(nbds$Name, getcount)

labs <- c(0,50,100,150)
pal.crime <- colorNumeric(colorRamp(c('#ccccff', 'red')), labs)

leaflet(nbds) %>% 
  addProviderTiles(providers$CartoDB.PositronNoLabels) %>%
  addPolygons(stroke=T,
              weight=1,
              popup=paste0(nbds$Name,"<br/>Shootings: ",nbds$count),
              color=~pal.crime(count),
              fillOpacity=.5) %>%
  addLegend("bottomright",title="# of Shootings (2012-present)",colors=~pal.crime(labs),labels=~labs)

Ceasefire weekends

Ceasefires have been called four times per year since August 2017. These are “ceasefire weekends” but their impact often extends well beyond a few days (one lasted twelve). Check out their website.

Here are the first days of each ceasefire (all Fridays).

# first day (Friday) of ceasefire weekends
ceasefire.initial <- as.Date(c("08/04/2017",
                               "11/03/2017",
                               "02/02/2018",
                               "05/11/2018",
                               "08/03/2018",
                               "11/02/2018",
                               "02/01/2019",
                               "05/10/2019",
                               "08/02/2019"),
                             format="%m/%d/%Y")

We determine the weekends of each ceasefire so we can see whether these weekends had fewer shootings, after accounting for other trends and seasonal effects.

ceasefire.weekends <- ceasefire.initial[1]
for (i in 2:length(ceasefire.initial)){
  ceasefire.weekends <- c(ceasefire.weekends,
                          seq(from=ceasefire.initial[i], 
                              by="day", 
                              length.out=3))
}

Initial model

First we’re going to use Facebook’s prophet package to check the rough effect of the ceasefire (coded as a “holiday”) while accounting for time trends and yearly and weekly seasonality. The prophet package was designed to be a good first pass: it gives you decent forecasts without a lot of manual effort.

We feed the ceasefire weekend dates into the model as special days or “holidays.”

ceasefires <- data_frame(
  holiday = 'ceasefire',
  ds = ceasefire.initial, # Fridays
  lower_window = 0,
  upper_window = 2 # for Sat and Sun
)

Then fit the model, accounting for general time trends as well as yearly and weekly seasonality.

library(prophet)

m1 <- prophet(data.frame(ds=daily$CrimeDate,
                         y=daily$shootings), 
              yearly.seasonality = T,
              weekly.seasonality = T, 
              mcmc.samples = 500,
              holidays = ceasefires,
              cores=4)

We plot the result. We can see that there are more shootings in recent years, but ceasefires seem to be associated with fewer shootings. 

future <- make_future_dataframe(m1, periods = 90)
forecast <- predict(m1, future)
plot(m1, forecast)

We can break the time trend down into components. The “holidays” are ceasefires, and suggest that the ceasefires decrease shootings even after accounting for weekly seasonality, yearly seasonality, and overall time trends. (Also, in general, it looks like Thursdays are good, and summers are bad.)

prophet_plot_components(m1, forecast)

Customized model

The prophet package doesn’t let us fit poisson likelihoods and doesn’t let us drill down and check the various effects for stuff like statistical significance, so we’re going to fit our own Bayesian model using Stan via the brms package.

We give the model information about the day of the week (Mon-Sun), the day of the year (0-364), and the date (julian), as well as a binary variable indicating whether a date occurs during ceasefire.

library(lubridate)

daily$weekday <- factor(weekdays(daily$CrimeDate), levels=c("Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday"))
daily$day.of.year <- yday(daily$CrimeDate)
daily$ceasefire <- factor(ifelse(daily$CrimeDate %in% ceasefire.weekends, 1, 0),labels=c("Regular Day","Ceasefire Weekend"))
daily$jul <- julian(daily$CrimeDate)
daily$julz <- scale(daily$jul)[,1]
Fit the model

We’ll predict shootings using a spline time trend, a cyclical spline for yearly seasonality, random effects for day of the week, and a binary indicator for the ceasefire. We use a Poisson link function because the outcome is a count and the mean is about the same as sd.

library(brms)

m2 <- brm(shootings ~ s(julz) +
            s(day.of.year, bs="cc") + #cyclical constraint 
             weekday + 
             ceasefire, 
           data=daily,
           cores=4,
           iter=500,
           family=poisson)
Plot model components

The results cohere well with the prophet model, but now we can drill down further. First the full model. The Ceasefires are visible as six dramatic downward spikes in 2017-2018.

And here are the decomposed seasonal and time trend effects.

daily$Estimate <- fitted(m2, scale="response")[,1]
# 80% posterior predictive interval
preds <- predict(m2, summary=FALSE)
daily$high <- apply(preds,2,function(x){quantile(x, prob=c(.9))})
daily$low <- apply(preds,2,function(x){quantile(x, prob=c(.1))})

daily %>% 
  ggplot(aes(x = CrimeDate, y = shootings)) +
  geom_point(alpha=.2) +
  geom_line(aes(y = Estimate), alpha=.5, color="red") +
  geom_ribbon(aes(ymin=low, ymax=high), alpha=.2) +
  xlab("date") +
  theme_bw()
p <- purrr::map(
       c("day.of.year [all]", "weekday", "julz [all]", "ceasefire"),
       ~  ggeffects::ggpredict(m2, .x) %>% plot() + theme_bw()
     )

p[[1]] <- p[[1]] + xlab("Day of year") + ggtitle(" ")
p[[2]] <- p[[2]] + xlab("Day of week") + ggtitle(" ")
p[[3]] <- p[[3]] + xlab("Time trend") + ggtitle(" ")
library(gridExtra)
grid.arrange(p[[3]], p[[2]],p[[1]])

And now look at the effect of ceasefire weekend, after accounting for these time trends and seasonal effects. On an average day in Baltimore, at least two people get shot. But ceasefire weekend is different. Ceasefire weekends appear to prevent more than one shooting every day.

p[[4]] + xlab("")